進階繪圖功能

林嶔 (Lin, Chin)

Lesson 10

第一節：3D圖形(1)

– 真的想要畫好看的圖，Google後修改前人的程式碼最快

F9_3

第一節：3D圖形(2)

F4_1

– 其數據集包含了150個樣本，都屬於鳶尾屬下的三個亞屬，分別是山鳶尾(setosa)、變色鳶尾(versicolor)和維吉尼亞鳶尾(virginica)。四個特徵被用作樣本的定量分析，它們分別是花萼(sepal)和花瓣(petal)的長度和寬度。

data(iris)
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

第一節：3D圖形(3)

– 我們先用2D散布圖來看看…

COLOR = as.integer(iris$Species)+1 #先根據Species分顏色，顏色代碼2在R裡面是紅色；3是綠色；4是藍色

par(mfrow = c(2, 3))
plot(iris[,"Sepal.Length"], iris[,"Sepal.Width"], pch = 19, col = COLOR)
legend("topright", levels(iris$Species), pch = 19, col = 2:4, bg = "gray90")
plot(iris[,"Sepal.Length"], iris[,"Petal.Length"], pch = 19, col = COLOR)
legend("bottomright", levels(iris$Species), pch = 19, col = 2:4, bg = "gray90")
plot(iris[,"Sepal.Length"], iris[,"Petal.Width"], pch = 19, col = COLOR)
legend("bottomright", levels(iris$Species), pch = 19, col = 2:4, bg = "gray90")
plot(iris[,"Sepal.Width"], iris[,"Petal.Length"], pch = 19, col = COLOR)
legend("topright", levels(iris$Species), pch = 19, col = 2:4, bg = "gray90")
plot(iris[,"Sepal.Width"], iris[,"Petal.Width"], pch = 19, col = COLOR)
legend("topright", levels(iris$Species), pch = 19, col = 2:4, bg = "gray90")
plot(iris[,"Petal.Length"], iris[,"Petal.Width"], pch = 19, col = COLOR)
legend("bottomright", levels(iris$Species), pch = 19, col = 2:4, bg = "gray90")

第一節：3D圖形(4)

library(scatterplot3d)
scatterplot3d(x = iris[,"Sepal.Length"],
              y = iris[,"Sepal.Width"],
              z = iris[,"Petal.Length"],
              color = COLOR, pch = 19, angle = 40, main="3D Scatterplot")

第一節：3D圖形(5)

– 套件『rgl』是在R裡面最常拿來繪製3D圖形的套件，他支援了互動式的3D圖像。

library(rgl)
library(rgl)
plot3d(x = iris[,"Sepal.Length"],
       y = iris[,"Sepal.Width"],
       z = iris[,"Petal.Length"],
       col = COLOR, size = 3, main="3D Scatterplot")

F9_4

練習1：學習如何做出進階功能

– 我們來牛刀小試一下，剛剛的圖片我們想要命令R將我們的點圈起來，我們可以使用下面程式碼…

library(rgl)
plot3d(x = iris[,"Sepal.Length"],
       y = iris[,"Sepal.Width"],
       z = iris[,"Petal.Length"],
       col = COLOR, size = 3, main="3D Scatterplot")

VCOV = cov(iris[,c("Sepal.Length", "Sepal.Width", "Petal.Length")])
MEAN = c(mean(iris[,"Sepal.Length"]), mean(iris[,"Sepal.Width"]), mean(iris[,"Petal.Length"]))

ellips = ellipse3d(VCOV, centre = MEAN, level = 0.95)
plot3d(ellips, col = "black", alpha = 0.2, add = TRUE, box = FALSE)

F9_5

You must enable Javascript to view this page properly.

– 如果你google的到，你甚至可以利用套件『rgl』做出gif的動畫檔案。(註：不一定要全程使用R做出來，只要能做出來就是好方法！)

F9_7

練習1答案

plot3d(x = iris[,"Sepal.Length"],
       y = iris[,"Sepal.Width"],
       z = iris[,"Petal.Length"],
       col = COLOR, size = 3, xlab = "Sepal Length", ylab = "Sepal Width", zlab = "Petal Length")

VCOV = cov(iris[1:50,c("Sepal.Length", "Sepal.Width", "Petal.Length")])
MEAN = apply(iris[1:50,c("Sepal.Length", "Sepal.Width", "Petal.Length")], 2, mean)

ellips <- ellipse3d(VCOV, centre = MEAN, level = 0.95)
plot3d(ellips, col = "red", alpha = 0.2, add = TRUE, box = FALSE)


VCOV = cov(iris[51:100,c("Sepal.Length", "Sepal.Width", "Petal.Length")])
MEAN = apply(iris[51:100,c("Sepal.Length", "Sepal.Width", "Petal.Length")], 2, mean)

ellips <- ellipse3d(VCOV, centre = MEAN, level = 0.95)
plot3d(ellips, col = "green", alpha = 0.2, add = TRUE, box = FALSE)


VCOV = cov(iris[101:150,c("Sepal.Length", "Sepal.Width", "Petal.Length")])
MEAN = apply(iris[101:150,c("Sepal.Length", "Sepal.Width", "Petal.Length")], 2, mean)

ellips <- ellipse3d(VCOV, centre = MEAN, level = 0.95)
plot3d(ellips, col = "blue", alpha = 0.2, add = TRUE, box = FALSE)

You must enable Javascript to view this page properly.

第二節：結合地理資訊的視覺化(1)

– Google地圖是我們最常用的地圖，他包含了街道圖/衛星圖等不同圖資，在R裡面我們可以透過套件輕鬆的與Google地圖作結合。

– 請至這裡下載登革熱流行的範例資料。

dat = read.csv("data10_1.csv", header = TRUE, fileEncoding = 'CP950')
head(dat)
##       發病日 個案研判日     通報日 性別 年齡層 居住縣市 居住鄉鎮 居住村里
## 1 1998/01/02            1998/01/07   男  40-44   屏東縣   屏東市         
## 2 1998/01/03            1998/01/14   男  30-34   屏東縣   東港鎮         
## 3 1998/01/13            1998/02/18   男  55-59   宜蘭縣   宜蘭市         
## 4 1998/01/15            1998/01/23   男  35-39   高雄市   苓雅區         
## 5 1998/01/20            1998/02/04   男  55-59   宜蘭縣   五結鄉         
## 6 1998/01/22            1998/02/19   男  20-24   桃園市   蘆竹區         
##      最小統計區 最小統計區中心點X 最小統計區中心點Y   一級統計區 二級統計區
## 1 A1320-0136-00          120.5059          22.46425 A1320-04-008   A1320-04
## 2 A1303-0150-00          120.4536          22.46639 A1303-09-007   A1303-09
## 3 A0201-0449-00          121.7514          24.74922 A0201-23-006   A0201-23
## 4 A6408-0153-00          120.3381          22.63032 A6408-10-010   A6408-10
## 5 A0209-0232-00          121.7983          24.68457 A0209-10-005   A0209-10
## 6 A0305-0543-00          121.2965          25.04431 A0305-36-002   A0305-36
##   感染縣市 感染鄉鎮 感染村里 是否境外移入 感染國家 確定病例數
## 1                                      否                   1
## 2                                      是                   1
## 3                                      是                   1
## 4                                      否                   1
## 5                                      否                   1
## 6                                      是                   1
install.packages("RgoogleMaps", repos="http://R-Forge.R-project.org")
library(RgoogleMaps)

第二節：結合地理資訊的視覺化(2)

– 注意，你需要準備API Key：

lat = c(22.88751, 23.41373)
lon = c(120.023, 120.6562)
center = c(mean(lat), mean(lon))
zoom = min(MaxZoom(range(lat), range(lon)))

MyMap = GetMap(center = center, zoom = zoom, API_console_key = 'AIzaSyA4DVFtF70aXE7RgrXViy2z5Ku2pMkVxFI')
PlotOnStaticMap(MyMap)

MyMap2 = GetMap(center = center, zoom = zoom, maptype = "satellite", API_console_key = 'AIzaSyA4DVFtF70aXE7RgrXViy2z5Ku2pMkVxFI')
PlotOnStaticMap(MyMap2)

第二節：結合地理資訊的視覺化(3)

dat[,1] = as.Date(dat[,1])
subdat = dat[dat[,1] <= as.Date("2015-09-30") & dat[,1] >= as.Date("2015-09-01") & dat[,6] == "台南市",]
nrow(subdat)
## [1] 12815

– 我們可以透過函數「PlotOnStaticMap()」把點放到MyMap上面

PlotOnStaticMap(MyMap, lat = subdat$最小統計區中心點Y, lon = subdat$最小統計區中心點X, pch = 19, col = "red", cex = 1)

第二節：結合地理資訊的視覺化(4)

x1 <- subdat$最小統計區中心點Y
x2 <- subdat$最小統計區中心點X
df <- data.frame(x1,x2)

## Use densCols() output to get density at each point
x <- densCols(x1,x2, colramp=colorRampPalette(c("black", "white")))
df$dens <- col2rgb(x)[1,] + 1L

## Map densities to colors
cols <-  colorRampPalette(c("#000099", "#00FEFF", "#45FE4F", 
                            "#FCFF00", "#FF9400", "#FF3100"))(256)
df$col <- cols[df$dens]

PlotOnStaticMap(MyMap, lat = df$x1, lon = df$x2, pch = 19, col = df$col, cex = 1.5)

第二節：結合地理資訊的視覺化(5)

mapdat = read.table("TWmap/TWmap/鄉村資料.txt", fileEncoding = "UTF-8")
head(mapdat)
##   objectid_1 village town county village_id  town_id county_id     tv_all
## 0          1  中山里 東區 嘉義市        061 10020010     10020 東區中山里
## 1          2  後驛里 西區 嘉義市        064 10020020     10020 西區後驛里
## 2          3  東川里 東區 嘉義市        038 10020010     10020 東區東川里
## 3          4  番社里 西區 嘉義市        061 10020020     10020 西區番社里
## 4          5  中央里 東區 嘉義市        062 10020010     10020 東區中央里
## 5          6  東興里 東區 嘉義市        060 10020010     10020 東區東興里
##       villcode ori_villag     area    max_x   max_y    min_x   min_y        x
## 0 A2001-061-00       NULL 230686.8 195185.2 2598108 194071.2 2597692 194698.8
## 1 A2002-064-00       NULL 363016.4 193222.2 2598067 192319.2 2597159 192755.7
## 2 A2001-038-00       NULL 630623.2 196096.4 2598009 194992.9 2596974 195507.5
## 3 A2002-061-00       NULL 170099.2 193546.0 2597889 192872.0 2597386 193249.6
## 4 A2001-062-00       NULL 194986.9 194170.0 2597852 193726.8 2597264 193954.0
## 5 A2001-060-00       NULL 215453.5 195207.0 2597837 194114.4 2597456 194621.2
##         y         v_id sort id countyname townname objectid orig_fid shape_leng
## 0 2597889 10020010-061   20  0       NULL     NULL        0        0          0
## 1 2597628 10020020-064   20  1       NULL     NULL        0        0          0
## 2 2597570 10020010-038   20  2       NULL     NULL        0        0          0
## 3 2597626 10020020-061   20  3       NULL     NULL        0        0          0
## 4 2597550 10020010-062   20  4       NULL     NULL        0        0          0
## 5 2597662 10020010-060   20  5       NULL     NULL        0        0          0
##   shape_le_1 shape_area
## 0   2734.726   230686.8
## 1   2617.877   363016.4
## 2   3891.870   630623.2
## 3   2021.745   170099.2
## 4   1812.928   194986.9
## 5   2658.314   215453.5
Tainan.mapdat = mapdat[mapdat[,4] == "臺南市",]
nrow(Tainan.mapdat)
## [1] 752
MysubMap = GetMap(center = center, zoom = zoom, maptype = "satellite", API_console_key = 'AIzaSyA4DVFtF70aXE7RgrXViy2z5Ku2pMkVxFI')
PlotOnStaticMap(MysubMap)

for (i in 1:nrow(Tainan.mapdat)) {
  linedat = read.csv(paste0("TWmap/TWmap/編號", Tainan.mapdat[i,1], ".csv"), header = TRUE, fileEncoding = 'CP950')
  PlotOnStaticMap(MysubMap, lat = linedat[,2], lon = linedat[,1], FUN = lines, add = TRUE, col = "red")
}

MysubMap = GetMap(center = center, zoom = zoom, maptype = "satellite", API_console_key = 'AIzaSyA4DVFtF70aXE7RgrXViy2z5Ku2pMkVxFI')
PlotOnStaticMap(MysubMap)

for (i in 1:nrow(Tainan.mapdat)) {
  linedat = read.csv(paste0("TWmap/TWmap/編號", Tainan.mapdat[i,1], ".csv"), header = TRUE, fileEncoding = 'CP950')
  PlotOnStaticMap(MysubMap, lat = linedat[,2], lon = linedat[,1], FUN = lines, add = TRUE, col = "red")
}

PlotOnStaticMap(MysubMap, lat = df$x1, lon = df$x2, pch = 19, col = df$col, add = TRUE)

第二節：結合地理資訊的視覺化(6)

lat = c(22.88751, 23.41373)
lon = c(120.023, 120.6562)
plot.new()
plot.window(xlim = lon, ylim = lat)

for (i in 1:nrow(Tainan.mapdat)) {
  linedat = read.csv(paste0("TWmap/TWmap/編號", Tainan.mapdat[i,1], ".csv"), header = TRUE, fileEncoding = 'CP950')
  polygon(linedat[,1], linedat[,2], col = "white")
}

lat = c(22.88751, 23.41373)
lon = c(120.023, 120.6562)
plot.new()
plot.window(xlim = lon, ylim = lat)

for (i in 1:nrow(Tainan.mapdat)) {
  linedat = read.csv(paste0("TWmap/TWmap/編號", Tainan.mapdat[i,1], ".csv"), header = TRUE, fileEncoding = 'CP950')
  n.sample = sum(subdat[,16] == as.character(Tainan.mapdat[i,2]))
  if (n.sample == 0) {COL = "#FFFFFF"} 
  else if (n.sample <= 3) {COL = "#000099"} 
  else if (n.sample <= 10) {COL = "#00FEFF"} 
  else if (n.sample <= 30) {COL = "#45FE4F"} 
  else if (n.sample <= 50) {COL = "#FCFF00"} 
  else if (n.sample <= 100) {COL = "#FF9400"} 
  else {COL = "#FF3100"} 
  polygon(linedat[,1], linedat[,2], col = COL)
}

第二節：結合地理資訊的視覺化(7)

lat = c(22.88751, 23.41373)
lon = c(120.023, 120.6562)
plot.new()
plot.window(xlim = lon, ylim = lat)

for (i in 1:nrow(Tainan.mapdat)) {
  linedat = read.csv(paste0("TWmap/TWmap/編號", Tainan.mapdat[i,1], ".csv"), header = TRUE, fileEncoding = 'CP950')
  n.sample = sum(subdat[,16] == as.character(Tainan.mapdat[i,2]))
  if (n.sample == 0) {COL = "#FFFFFF80"} 
  else if (n.sample <= 3) {COL = "#00009980"} 
  else if (n.sample <= 10) {COL = "#00FEFF80"} 
  else if (n.sample <= 30) {COL = "#45FE4F80"} 
  else if (n.sample <= 50) {COL = "#FCFF0080"} 
  else if (n.sample <= 100) {COL = "#FF940080"} 
  else {COL = "#FF310080"} 
  polygon(linedat[,1], linedat[,2], col = COL)
}

MysubMap = GetMap(center = center, zoom = zoom, maptype = "satellite", API_console_key = 'AIzaSyA4DVFtF70aXE7RgrXViy2z5Ku2pMkVxFI')
PlotOnStaticMap(MysubMap)

for (i in 1:nrow(Tainan.mapdat)) {
  linedat = read.csv(paste0("TWmap/TWmap/編號", Tainan.mapdat[i,1], ".csv"), header = TRUE, fileEncoding = 'CP950')
  n.sample = sum(subdat[,16] == as.character(Tainan.mapdat[i,2]))
  if (n.sample == 0) {COL = "#FFFFFF80"} 
  else if (n.sample <= 3) {COL = "#00009980"} 
  else if (n.sample <= 10) {COL = "#00FEFF80"} 
  else if (n.sample <= 30) {COL = "#45FE4F80"} 
  else if (n.sample <= 50) {COL = "#FCFF0080"} 
  else if (n.sample <= 100) {COL = "#FF940080"} 
  else {COL = "#FF310080"} 
  PlotOnStaticMap(MysubMap, lat = linedat[,2], lon = linedat[,1], FUN = polygon, add = TRUE, col = COL)
}

練習2：做出逐日變化的動畫

– 假定你希望呈現整個9月份『逐日』的『累積個案散布圖』，你有辦法熟練得操作你的畫圖技術了嗎?

F10_1

練習2答案

MysubMap = GetMap(center = center, zoom = zoom, maptype = "satellite", API_console_key = 'AIzaSyA4DVFtF70aXE7RgrXViy2z5Ku2pMkVxFI')

pdf('Tainan.pdf', width = 7, height = 7)

for (k in 1:30) {
  subdat = dat[dat[,1] <= as.Date(paste0("2015-09-", k)) & dat[,1] >= as.Date("2015-09-01") & dat[,6] == "台南市",]
  PlotOnStaticMap(MysubMap)
  for (i in 1:nrow(Tainan.mapdat)) {
    linedat = read.csv(paste0("TWmap/TWmap/編號", Tainan.mapdat[i,1], ".csv"), header = TRUE, fileEncoding = 'CP950')
    n.sample = sum(subdat[,16] == as.character(Tainan.mapdat[i,2]))
    if (n.sample == 0) {COL = "#FFFFFF80"} 
    else if (n.sample <= 3) {COL = "#00009980"} 
    else if (n.sample <= 10) {COL = "#00FEFF80"} 
    else if (n.sample <= 30) {COL = "#45FE4F80"} 
    else if (n.sample <= 50) {COL = "#FCFF0080"} 
    else if (n.sample <= 100) {COL = "#FF940080"} 
    else {COL = "#FF310080"} 
    PlotOnStaticMap(MysubMap, lat = linedat[,2], lon = linedat[,1], FUN = polygon, add = TRUE, col = COL)
  }
  text(-200, -250, as.Date(paste0("2015-09-", k)), cex = 2, col = "red")
}

dev.off()

小結

– 透過一系列的學習，加上網路上的資料，資料整理、分析、呈現我們都已經有初步了解，下一步就是資料蒐集的部分，我們有沒有可能也自動化的蒐集資料呢?